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Drexel  University 
3141  Chestnut  St. 

Philadelphia,  PA  19104 

ABSTRACT 

Rational  design  and  selection  of  candidate  porous  biomaterials  to  serve  as  tissue 
engineering  constructs  rests  on  our  ability  to  understand  the  influence  of  the  porous 
microarchitecture  on  the  transport  of  chemical  species  (e.g.,  nutrients  and  signaling 
compounds),  fluid  flow,  and  cellular  locomotion  and  growth.  We  have  begun  to  study  the 
behavior  of  chemotactically  mobile  cells  in  response  to  unsteady  signaling  molecule 
concentration  fields  using  a computational  simulation-based  model.  The  model  couples 
fully  time-dependent  finite-difference  solution  of  a reaction-diffusion  equation  for  the 
concentration  field  of  a generic  chemoattractant,  to  biased  random  walks  representing 
individual  moving  cells.  This  model  is  a first  step  in  building  a quantitative,  pore-level 
model  of  mass  and  cellular  transport  in  porous  tissue-engineered  constructs.  In  these 
proceedings,  we  focus  on  our  recent  findings  regarding  the  influence  of  flux-reactive 
boundary  conditions  in  heterogeneous  2D  domains  on  the  chemotactic  response  of  otherwise 
randomly  moving  cells.  In  particular,  we  find  that,  when  cells  are  forced  to  “crawl”  around 
obstacles  in  order  to  approach  a point  source  of  chemoattractant,  the  reactivity  of  the 
obstacle  surface  with  respect  to  the  chemoattractant  strongly  determines  the  morphology 
of  the  cells’  paths  of  locomotion.  Cells  crawl  along  non-reactive  surfaces  and  strongly  avoid 
reactive  surfaces,  due  to  the  nature  of  the  chemoattractant  concentration  gradients  near 
the  surface.  We  show  further  that  tuning  the  reactivity  of  the  surfaces  of  two  obstacles 
defining  a gap  can  control  the  passage  of  cells  through  the  gap.  From  our  work,  we  infer 
the  importance  of  a proper  treatment  of  boundary  conditions  in  any  future  pore-level 
quantitatve  modeling  of  mass  transport  and  cellular  response  in  porous  media. 

INTRODUCTION 

Current  technology  gives  us  the  ability  to  construct  a wide  variety  of  porous 
biomaterials  for  therapeutic  applications,  including  tissue  engineering  and  immuno-isolated 
implants.  In  fact,  this  ability  far  outpaces  our  understanding,  at  a fundamental  level,  of 
why  one  porous  material  performs  better  than  another  in  any  given  application.  In  many 
applications,  one  of  the  most  important  performance  requirements  is  that  the  material 
allow  for  easy  development  of  healthy,  internal  capillary  networks  [1] . This 
neovascularization  within  an  implant  greatly  reduces  the  impediment  to  biochemical 
transport  relative  to  diffusion  through  the  otherwise  present  fibrous  capsule  characteristic 
of  the  classical  foreign  body  response  (FBR)  [2].  Mitigating  the  FBR  in  this  way  allows  for 
optimal  nutrient /waste  transport  for  engineered  tissues  grown  on  porous  scaffolds  [3]  and 
has  the  potential  to  guarantee  rapid,  controllable  dosing  in  drug  delivery  and 
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Figure  1.  Schematic  of  capillary  sprout  growth  in  angiogenesis,  showing  a close-up  of  the 
tip  in  longitudinal  cross-section. 


iinmuno-isolated  gene  and  cell-based  implant  therapies  [4],  It  is  therefore  vitally  important 
that  we  explore  ways  to  control  this  neovasculature  development.  One  goal  must  be 
incorporating  this  ability,  which  we  refer  to  as  “directed  angiogenesis,”  into  design 
strategies  for  porous  biomaterials.  This  requires  in  part  a fundamental  understanding  of 
the  interactions  between  porous  microstructures  and  growing  capillary  sprouts. 

The  directed  migration  of  endothelial  cells  (ECs)  is  the  mechanism  underlying  the 
construction  of  capillary  networks  as  a response  to  stimulants  such  as  vascular  endothelial 
growth  factor  (VEGF)  and  tumor  angiogenic  growth  factor  (TAF)  [5].  The  morphology  of 
a sprout  is  determined  by  (i)  the  biased  locomotion  path  of  the  single  EC  at  the  sprout  tip, 
and  (ii)  the  proliferation  of  cells  in  the  tubule  following  the  sprout  tip,  ns  depicted  in 
figure  1.  The  EC  locomotion  is  directed  over  large  length  scales,  while  at  small  length 
scales  it  is  apparently  random,  owing  to  inhomogeneities  in  the  ECM  and  void  space 
through  which  these  cells  must  crawl,  as  well  as  concentration  variations  in  stimulants  and 
nutrients,  among  others.  This  biased  locomotion  is  categorized  as  a mixture  of 
“chemotaxis,”  in  which  the  direction  of  cell  motion  is  influenced  by  external  concentration 
gradients,  and  “chemokinesis,”  in  which  the  speed  of  cell  motion  is  influenced  by  these 
gradients. 

There  is  evidence  in  the  experimental  literature  that,  chemotaxis  can  be  harnessed  to 
induce  EC  penetration  of  porous  scaffolds  which  would  normally  not  allow  for  EC 
penetration.  Various  biocompatible  porous  scaffold  materials,  including  expanded 
poly(tetrafluoroethylene)  (ePTFE)  [6]  and  poly  (vinyl  alcohol)  [3],  have  been  shown  to 
admit  EC  penetration  in  vitro  when  the  mean  pore  size  is  greater  than  ~60  pm.  However, 
Kidd  et  al.  showed  that  when  an  ePTFE  disc  is  treated  with  cellular  detritus  rich  in 
growth  factors  on  one  side,  microvcsscl  penetration  is  possible  even  through  pores  as  small 
as  30  pm  [7].  A similar  pretreatment  strategy  was  employed  by  Sanders  et  al.,  who 
demonstrated  the  feasibility  of  “prevascularizing”  a polyurethane  mesh  biomaterial  implant, 
by  placing  it  in  temporary  contact  with  the  relatively  active  vasculature  of  the 
chorioallontoic  membrane  of  quail  embryos  [8],  The  quail  vasculature  readily  penetrated  the 
implant,  which  was  subsequently  removed  from  contact  with  the  membrane  and  implanted 
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in  rats.  This  means  that  both  sets  of  authors  inferred  that  concentration  gradients  of  some 
type  of  signaling  molecule  influenced  the  direction  and  location  of  vessel  growth. 

The  great  interest  in  the  biomaterials  community  in  directed  angiogenesis  has  led  to  a 
variety  of  empirical  optimization  approaches.  The  biomaterials  engineer  can  imagine  and 
construct  a variety  of  materials  with  different  microporous  architectures,  can  biochemically 
modify  these  structures,  and  design  any  number  of  protocols  for  ex  vivo  culturing.  In  the 
face  of  this  variety,  we  must,  place  rational  limits  on  the  amount  of  trial-and-error,  and 
therefore  the  expense,  required  to  decide  which  porous  biomaterial  design  strategy  best 
suits  a particular  application.  A significant  step  forward  in  this  regard  would  be  a 
simulation  model  which  allows  for  rapid  in  virtuo  testing  of  a proposed  microarchitecture, 
together  with  a proposed  biochemical  pretreatment  and  culturing  protocol,  prior  to  its 
fabrication  and  in  vitro/in  vivo  testing.  Such  a model  must,  at  a minimum,  capture  the 
3-D  microstructure  of  a porous  material  and  account  for  chemotaxis  and  chemokinesis  in 
cell  migration,  to  quantify  effects  of  length  scales  and  correlations  in  the  porous  structure 
and  biochemical  pretreatment  techniques  on  the  rate  of  cell  penetration. 

No  such  model  yet  exists.  However,  our  work  toward  such  a model  has  led  us  to 
investigate  chemotactic  locomotion  in  well-controlled  simulations  in  inhomogeneous  2D 
domains.  Living  cells  are  too  complicated  to  model  with  complete  accuracy.  In 
multicellular  simulations,  simple  rules  governing  cell  behavior  are  imposed,  and  emergent 
behavior  from  cell  populations  is  observed  and  analyzed.  A “cell”  is  a computational  entity 
that  can  move,  divide,  die,  perhaps  transform  itself,  consume  nutrients,  and  produce 
wastes.  The  task  of  the  simulation  developer  is  to  come  up  with  sensible  rules  governing 
this  behavior,  so  that,  to  an  observer,  the  cell  mimics  the  behavior  of  a living  cell.  This 
non-trivial  task  means  that  several  dozen  parameter  values  must  be  carefully  assigned 
based  on  existing  experimental  results,  and  in  some  instances,  free  parameters  are  fit  to 
experimental  results  in  the  process  of  building  the  model. 

A popular  method  in  multicellular  simulations  is  cellular  automata  (CA).  In  CA,  an 
automaton  is  an  array  of  “cells”  (not  in  the  biological  sense) , each  being  a box  in  space 
that  can  have  one  of  many  states.  In  biological  applications,  a cell  can  be  unoccupied,  or 
occupied  with  a computational  cell.  In  each  simulation  pass  or  time-step,  each  occupied 
cell  is  given  the  chance  to  change  its  state  based  on  its  local  environment.  Relatively 
simple  rules  for  division  and  death  result  in  surprisingly  rich  and  structured  behaviors,  as 
anyone  who  has  watched  J.  H.  Conway’s  “Game  of  Life”  screensaver  can  attest  [9].  A 
relatively  more  complex  set  of  rules  was  used  to  model  growing  endothelial  cells  in 
pioneering  work  by  Zygourakis  et  al.  (e.g.,  [10,  11]).  Quite  recent  work  has  used  CA  in 
conjunction  with  simple  transport  models  to  probe  the  effects  cell-to-cell  signaling  in  the 
development  of  microvascular  structures  [12]. 

A related  technique  in  multicellular  simulations  is  the  biased  random  walk 
(BRW)  [13,  14,  15].  This  type  of  model  has  found  widespread  use  in  the  simulation  of 
angiogenesis.  This  path  may  generally  tend  toward  a source  of  chemoattractant,  but  also 
has  some  tunable  random  contribution.  Most  studies  to  date  have  assumed  that  the 
concentration  field  of  chemoattractant  is  steady  in  time.  The  only  true  distinction  between 
CA  and  BRW  cell  models  is  that.  CA  restricts  cell  positions  to  lattice  sites,  which  BRW  are 
typically  off-lattice.  In  only  one  instance  has  a BRW  multicellular  simulation  been  coupled 
with  the  fully  transient  solution  of  a reaction-diffusion  equation  for  chemoattractant  [14], 
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and  this  was  done  in  a 2D  homogeneous  domain.  This  is  a significant  advance  because  it 
dynamically  coupled  the  concentration  of  chemoattractant  and  migrating  cells.  Tong  and 
Yuan  [14]  demonstrated  that  this  interplay  is  responsible  for  the  well-known 
“brush-border’'  effect,  in  which  a very  high  density  of  capillary  branching  occurs  near  the 
source  of  a growth  factor,  screening  the  parent  vessel  from  the  growth  factor  and 
ultimately  slowing  the  total  process  of  angiogenesis. 

We  have  developed  a BRW  multicellular  simulation  in  the  same  spirit  of  Tong  and 
Yuan,  with  an  eye  toward  the  eventual  development  of  a multicelluar/transport  simulation 
model  of  tissue  ingrowth  in  porous  constructs.  The  purpose  of  this  contribution  is  to 
discuss  our  recent  work  investigating  the  effect  of  surfaces  winch  consume  chemoattractant 
on  the  chemotactic  response  of  crawling  cells  in  heterogeneous  2D  domains.  The  goal 
generally  is  to  understand  how  both  the  geometry  and  reactivity  of  surfaces  in  a structured 
domain  modulate  concentration  fields  of  stimulants  (chemoattractants)  and  how  these 
modulations  can  be  engineered  to  control  cell  locomotion. 

MODEL  AND  SIMULATIONS 


Our  model,  described  in  detail  elsewhere  [16],  includes  (1)  a fully  time-resolved  finite 
difference  solution  winch  governs  the  transport  of  chemoattractant,  and  (2)  a “biased” 
random  walk  simulation  of  mobile  cells  which  couples  the  direction  of  cell  motion  to  the 
simultaneously  evolving  concentration  field  of  chemoattractant.  The  model  considers  the 
following  fundamental  events:  (i)  release  of  a chemoattractant  from  sources  or  cells;  (ii) 
diffusion  of  chemoattractant  in  the  domain;  (iii)  randomness  in  the  direction  of  cell 
migration,  presumably  due  to  other  underlying  variables  which  yet  we  do  not  have  enough 
information;  (iv)  chemotactic  response  of  cells. 

TVansport  of  chemoattractant 


The  major  task  of  the  finite  difference  solver  is  to  obtain  the  concentration  of 
chemoattractant  as  a function  of  position  and  time,  c(r,  t)  (boldface  denotes  a vector 
quantity).  The  governing  transport  equation  is  as  follow's: 


dc 

dt 


dv2c + y 


<5(r  - r,), 


(1) 


where  ® is  the  diffusivity  of  chemoattractant  and  r,  is  the  location  of  the  i’th  point  source. 
Point  sources  produce  at  a rate  P(f).  vc  is  the  volume  of  a cell,  so  represents  the  rate 
at  which  the  concentration  of  chemoattractant  is  introduced  by  a point  source  at.  position 
r,.  S(x)  is  the  Dirac  delta  function: 


x = 0 
x^O 


(2) 


For  simplicity  at  this  stage  of  model  development,  cell  division,  uptake  and  degradation  of 
chemical  species  are  neglected.  Numerical  values  of  the  relative  quantities  in  this  model 
appear  in  Table  I.  Eq.  1 is  solved  numerically  using  a finite-difference  technique  referred  to 


40 


as  the  “hopscotch  method”  [17,  18].  The  2D  homogeneous  solution  domain  measures 
2 mm  x 2 mm,  as  shown  in  Fig.  2,  and  is  discretized  onto  a mesh  of  200x200  nodes.  The 
origin  of  the  domain  is  the  lower-left  corner,  and  the  center  point  of  the  domain  is  node 
(100,100).  From  here  on,  lengths  are  reported  in  units  of  mm  unless  otherwise  stated. 

The  following  initial  and  boundary  conditions  on  the  solution  c(r,  t ) are  enforced: 

c(r,0)  = 0 (3) 

— DVc  = kc  (4) 


Eq.  4 represents  reactive  boundary  conditions  applied  along  boundary  curves  in  the 
domain.  The  proportionality  constant  k measures  the  reactivity  of  the  boundaries.  For  this 
early  stage  study,  we  have  assumed  first  order  kinetics.  As  is  standard  in  these  types  of 
“reaction-diffusion”  problems,  the  dimensions  of  the  constant  k are  length/time;  we  choose 
to  measure  k in  units  of  cm/s.  The  various  boundary  curve  geometries  considered  (defining 
the  heterogeneous  domains)  is  described  in  the  “Simulation  Protocols”  section. 

The  dimensionless  parameter  which  measures  whether  diffusion  or  reaction  dominates 
the  boundary  condition  given  in  Eq.  4 is  the  Peclet  number: 


Pe  = 


kh 
D ’ 


(5) 


where  we  have  chosen  the  lattice  spacing  h as  the  relavent  length  scale.  For  the  parameters 
given  in  table  I (D  = 10~6  cm2/s,  h = 10~3  cm),  in  addition  to  considering  the  case  of 
k = 0 (i.e.,  “no-flux”  boundaries),  we  have  chosen  to  sweep  values  of  k from  a “low”  Pe  of 
1CT3  to  a “high”  Pe  of  100. 

As  presented,  the  only  sink  for  chemoattractant  so  far  discussed  is  reaction  at  boundary 
curves,  controlled  by  the  reactivity  k.  Additionally,  we  employ  “global”  boundary 
conditions  of  Dirichlet  type;  these  too  are  discussed  in  the  “Simulation  Protocols” 
section. 


Biased  random  walk  simulation  of  chemotaxis 


Random  walks  are  propagated  in  the  following  way.  We  first  stipulate  that  a walker 
make  a step  on  average  with  a frequency  /p,  and  that  step  lengths  are  randomly  chosen 
from  a Gaussian  distribution  with  mean  (A l)  and  standard  deviation  ctai-  fr  and  (A/)  are 
constrained  relative  to  one  another  to  guarantee  that,  over  long  times,  cells  move  with  an 
average  velocity,  Vq,  observed  experimentally.  We  freely  choose  fT  and  <7a(,  but  these 
values  (shown  in  table  I)  are  the  same  for  all  simulations  discussed  here.  For  each  mesh 
update  (every  2 s of  integration  time)  each  cell  in  the  simulation  is  given  the  opportunity 
to  move.  The  probability  that  a cell  will  take  a step  is  the  probability  of  observing  an 
event  which  occurs  with  a frequency  fT  in  a time  interval  At:  PT  = /pAf.  If  a uniform 
random  variate  between  0 and  1 is  chosen  less  than  Pp,  the  cell  under  consideration  takes  a 
step.  The  direction  of  a step,  denoted  by  the  unit  vector  p,  is  computed  as 


r(f)  - r (t  - Afceii)  Vc 
“|  r(f)-r(f-  Atceii)  | HVc| 


+ (1  — a — ft)  A. 


(6) 
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Table  I.  Nomenclature;  Baseline  values  or  ranges  of  important  model  parameters. 


Parameter 

AfCone 

h 

T> 

P 

fr 


Vmax 

(AZ) 


cr&i 

a 

3 


rc 

Vc 

k 


Description 

Time  step  for  mesh  update 
Mesli  discretization 
Growth  factor  diffusivity 
Chemoattractant  production  rate  per  cell 
Average  frequency  of  random  walk  updates 
Maximum  cell  speed 
Mean  step  length 

Width  of  Gaussian  step  length  distil. 
Step  persistence  parameter 
Step  chemotactic  parameter 
Cell  radius 
Cell  volume 

Rate  constant  for  boundary  reactivity 


Value  or  Range 
2 s 

10  /mi 

10-G  cm2/s  [14] 

1 Pg/h  [19] 
0.0167  s'1 
20  pm/li  [13] 
0.333  pin 
0.083  pm 
0 

0-1 

10  pm  [14] 

« Jjrrj! 

0 and  10  fi-0.1  cm/s 


The  vector  p is  the  weighted  sum  of  three  unit  vectors.  The  first  term,  which  is  weighted 
by  a,  defines  persistence,  representing  the  fact  that  mobile  cells  keep  a “memory”  of  their 
previous  step  direction.  The  second  term,  which  is  weighted  by  3,  represents  the 
chemotactic  response.  The  final  term,  carrying  the  remainder  of  the  weight,  is  a random 
direction,  A,  chosen  uniformly.  We  interpret  this  “random”  component  of  the  direction  as  a 
lumping  together  of  all  underlying  processes  governing  cell  migration  for  which  we  have  no 
information  or  purposely  ignore  for  simplicity.  Also  for  reasons  of  simplicity  in  this  early 
study,  we  have  chosen  to  neglect  persistence,  setting  a to  0 for  all  cases. 

The  motion  of  cells  are  coupled  to  the  transport  of  chemoattractant  because  a cell  must 
“measure”  a local  concentration  gradient.  Because  cell  positions  are  off-lattice,  a cell  at  an 
arbitrary  location  in  2D  space  samples  c and  Vc  by  linear  interpolation  among  the  four 
nodes  nearest  it.  At  each  cell  position  update,  the  instant  values  of  average  c and  Vc  are 
used  in  the  determination  of  p.  When  a cell  injects  chemoattractant  into  the  domain,  it 
likewise  partitions  the  amount  over  the  four  nodes  nearest  it  via  linear  interpolation. 

An  important  component  of  the  model  is  the  pair  of  threshold  conditions  required  for 
chemotactic  response.  In  optimal  conditions,  cells  can  detect  a 1 percent  difference  in 
concentration  of  chemoattractant  along  their  length  [20].  It  has  also  been  shown  in 
experiments  that  high  chemoattractant  concentration  results  in  cell  movements  which  are 
completely'  decorrelated  from  the  direction  of  a chemoattractant,  gradient  [21].  Therefore, 
we  prescribe  that  the  coefficient  3 also  serve  as  a switch  on  the  chemotactic  response.  For 
any  particular  cell's  average  c and  Vc,  the  value  of  3 used  to  determine  the  direction  of  its 
next  step  is  set  to  0 if  ^,crc  < 0.005,  which  we  refer  to  as  the  “gradient  threshold.”  We  also 
employ  a “concentration  threshold."  setting  3 = 0 when  c.  < 10~10  g/cm3  [19]. 
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Figure  2.  Schematic  simulation  domains  for  the  four  types  of  chemotaxis  simulations  con- 
sidered here,  (a)  Single  walkers  around  single  obstacle,  (b)  Dual  walkers  around  single 
obstacles,  (c)  Dual  walkers  around  dual  obstacles. 


Simulation  Protocols 

Single  walks  around  single  obstacles 

We  first  consider  the  simple  scenario  in  which  an  immobile  point  source  at  the  center 
point  produces  chemoattractant  at  a constant  rate  (P  = 1 pg/h).  Cells  are  launched  from 
a distance  of  0.5  mm  from  the  source  for  randomly  chosen  values  of  polar  angle.  6.  A 
circular  “obstacle”  of  diameter  0.2  mm  is  located  0.3  mm  to  the  left,  of  the  source 
(figure  2a).  A global  Dirichlet-type  boundary  condition,  c(R,t)  = 0,  is  enforced  at  R = 

1 mm.  100,000  walks  are  conducted  for  various  values  of  chemotacticity  /3  and  surface 
reactivity  k.  The  primary  observable  is  the  time  required  for  the  cells  to  traverse  the 
0.5  mm  distance  from  the  launching  radius  to  the  source,  referred  to  as  the 
“time-to-contact,”  tc,  as  a function  of  angular  position,  0. 

Dual  walks  around  single  obstacles 

We  next  consider  the  scenario  in  which  two  cells  are  placed  initially  0.4  mm  apart  with 
the  midpoint  of  the  segment  defined  by  their  positions  lies  on  the  center  point.  Centered 
on  the  center  point  is  a circular  obstacle  of  diameter  0.2  mm  (figure  2b).  The  same  global 
Dirichlet  boundary  condition  from  the  previous  case  is  employed.  The  two  cells  produce 
chemoattractant  in  pulses  defined  by  a duty  cycle:  secretion  at  a rate  of  6 P for  a random 
“on”  time,  t0 n,  selected  uniformly  in  the  range  10  min  < £on  < 30  min,  followed  by  no 
secretion  for  a random  “off”  time,  £0ff,  selected  uniformly  in  the  range  50  min  < f0ff  < 

150  min.  (We  demonstrated  previously  that  this  random  “pulsed”  secretion  of 
chemoattractant  guarantees  that  two  cells  can  successfully  signal  one  another  in  a 
homogeneous  domain  [16].)  100  statistically  equivalent  individual  dual-cell  simulations  are 
performed  for  various  values  of  j3  and  k.  The  primary  observable  is  the  time  required  for 
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Figure  3.  (a)  Distributions  in  time-to-contact  over  0 , for  single  walkers  in  the  domain  with 
the  single  obstacle,  for  various  obstacle  boundary  reactivities  k;  units  of  k in  the  figure  are 
10" 3 cm/s.  (b)  Mean  of  each  distribution  in  (a),  (c)  Example  locomotion  paths  for  selected 
reactivities;  0 = 0.30. 


the  two  cells  to  reach  one  another,  again  referred  to  as  the  “time-to-contact,”  tc. 

Dual  walks  around  dual  obstacles 

The  next  scenario  we  consider  is  similar  to  the  previous,  except  that  we  now  place  two 
circular  obstacles  of  diameter  0.2  mm  separated  by  a distance  D,  with  their  midpoint  on 
the  center  point.  The  quantity  L = D - 0.2  nun  is  referred  to  as  the  “gap  spacing.”  as  this 
is  the  minimal  distance  between  the  two  bounding  curves  of  the  obstacles.  Dual  cells  are 
launched  at  180°  opposition  such  that,  in  order  to  reach  one  another,  they  must  crawl 
through  the  gap  between  the  obstacles  (figure  2c).  We  examined  the  effects  of  both  surface 
reactivity  k and  gap  spacing  L on  the;  chemotaet.ie  signaling  and  response  of  the  two  cells. 
Here,  the  primary  observable  is  whether  particular  values  of  k and  L permit,  successful 
passage  by  one  or  both  cells.  Our  aim  is  to  make  a first  attempt  at  understanding  the 
synergistic  roles  of  both  pore  size  and  surface  reactivity  on  chemotaet.ie  response  in  a 
simple  model  system. 

RESULTS  AND  DISCUSSION 

Single  walks  around  single  obstacles 

In  figure  3a,  we  present  the  time-to-contact  distributions  as  a function  of  angle  0 for  the 
case  of  a single  obstacle  offset  from  the  origin.  Each  curve  represents  the  results  of  runs 
with  different  obstacle  boundary  reactivity,  k.  First,  note  the  expected  result  that,  in 
general,  for  cells  launched  in  the  vicinity  of  9 = n,  for  which  the  point  source  at  the  origin 
is  eclipsed  by  the  obstacle,  the  contact  time  is  higher  than  the  average.  We  see  that,  the 
more  reactive  the  surface,  the  more  influence  is  felt  by  the  contact  time  distribution.  We 
show  the  mean  of  each  distribution  in  figure  3b.  The  contact  time  appears  most  sensitive 
to  k in  the  range  10~5  cm/s  < k < 10-3  cm/s.  The  saturation  at  high  reactivities  is 
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Figure  4.  (a)  Mean  contact  time  for  two  cells  forced  to  move  around  an  obstacle  vs. 
obstacle  boundary  reactivity,  k,  to  secreted  chemoattractant,  for  various  chemotacticities,  /?. 
(b)  Example  locomotion  paths  from  the  dual-cell,  single-obstacle  runs,  for  various  obstacle 
boundary  reactivities:  0 = 0.30.;  values  of  k in  the  figure  are  in  units  of  10-3  cm/s 


unambiguously  determined  by  the  global  boundary  condition  at  r = R.  As  can  be  seen  in 
figure  3c,  for  highly  reactive  surfaces,  cells  initially  move  away  from  both  the  obstacle  and 
the  point  source.  This  is  because  Vc  always  points  away  from  reactive  boundaries.  Because 
of  the  global  boundary,  there  must  exist  a radius  at  which  Vc.  in  the  r direction  in  the 
vicinity  of  6 = n changes  sign.  In  this  region,  the  gradient  threshhold  is  violated  and  cells 
“lose  signal”  for  some  time  until  they  can  diffuse  into  a region  in  which  some  gradients  are 
detectable.  These  gradients  invariably  lead  the  walkers  around  the  obstacle  to  the  source. 
We  learn  from  this  small  study  that,  because  reactive  surfaces  cause  gradients  in  c which 
point  away  from  the  surface,  chemotactically  active  cells  tend  to  avoid  reactive  surfaces, 
even  if  that  leads  them  further  from  a point  source  of  chemoattractant. 

Dual  walks  around  single  obstacles 

In  figure  4a,  we  show  the  mean  contact  times  for  two  chemotactically  active  mobile 
cells,  each  producing  chemoattractant  using  the  pulsed  protocol  described  previously,  vs. 
obstacle  boundary  reactivity  for  various  values  of  chemotacticity.  For  each  k and  0,  100 
individual  simulations  contribute  to  each  average.  We  observe  a similar  trend  as  that  seen 
for  the  single  walkers,  namely,  that  the  greatest  sensitivity  to  k occurs  in  the  range 
10“5  cm/s  < k < 10-3  cm/s.  Furthermore,  the  stronger  the  chemotacticity,  the  smaller  is 
the  influence  of  the  boundary  reactivity  on  the  mean  contact  time.  We  also  show  example 
locomotion  paths  in  figure  4.  This  illustrates  the  reason  for  this  sensitivity:  for  large 
reactivities,  cells  initially  move  away  from  the  obstacle  and  each  other,  until  feeling  the 
influence  of  the  global  boundary,  at  which  time  they  are  led  to  one  another  while 
“confined”  between  the  obstacle  at  the  origin  and  the  global  boundary.  The  stronger  the 
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Figure  5.  (a)  “Phase-diagram”  of  successful  passage  of  signaling,  ehemotaetically  responsive 
cells  through  a gap  of  size  L defined  by  two  circular  obstacles  with  boundaries  with  reactivity 
k.  Solid  points  denote  the  region  in  (L.fc)-spnee  for  which  passage  is  allowed,  and  open 
symbols  those  regions  where  passage  is  not  allowed.  In  all  cases,  0 = 0.30.  (b)  Example 
locomotion  paths  for  a gap  spacing  of  0.4  mm  and  k = 0.  Horizontal  line  segments  denote 
the  launch  positions,  (c)  Same  as  (b),  with  k = 1.0. 


obstacle  boundary  reactivity,  the  closer  are  the  cells  pushed  toward  the  global  boundary 
while  moving  toward  one  another.  These  results  again  demonstrate  that  ehemotaetically 
active  mobile  cells  respond  predictably  to  the  presence  of  obstacles  with  surfaces  that 
consume  chemoattractant. 

Dual  walks  around  dual  obstacles 

We  now  turn  to  our  simulations  of  two  mobile  ehemotaetically  active  cells  separated  by 
a gap  defined  by  two  circular  obstacles  (figure  2c).  The  two  cells  on  either  side  of  the  gap 
secrete  chemoattractant  in  pulses  describing  a random  duty  cycle  which  is  known  to  result 
in  successful  ehemotaxis  in  a homogeneous  domain.  Our  objective  was  to  determine  how 
both  the  gap  spacing  L and  the  surface  reactivity,  k,  both  dictate  whether  one  or  both  cells 
can  enter  or  pass  through  the  gap.  In  figure  5a.  we  report  in  phase-diagram  format  the 
results  of  our  dual-walker  simulations  at  various  gap  spacings  L and  obstacle  boundary 
reactivities,  k.  Solid  points  represent  those  values  of  L and  k for  which  one  or  both  of  the 
cells  pass  through  the  gap  sucessfully,  while  the  open  symbols  refer  to  those  cases  for  which 
no  successful  passage  through  the  gap  is  observed.  We  have  considered  gap  spacings 
generally  between  0.2  and  1.2  mm,  much  larger  than  an  EC  (~  10  pm)  or  the  minimal  pore 
size  through  which  ECs  are  known  to  pass  (~  60  pm).  For  low  k,  no  hindrance  to  passage 
is  observed  for  any  gap  size  considered;  an  example  of  locomotion  paths  for  k = 0 and  L = 
0.4  mm  is  shown  in  figure  5b.  However,  at  moderate  to  high  reactivities,  gap  sizes  of  0.2  - 
0.5  mm  prevent  passage;  an  example  of  locomotion  paths  for  for  k = 10~3  cm/s  and  L = 
0.4  mm  is  shown  in  figure  5c.  One  can  see  that  the  cells  arc  simply  forced  away  from  the 
gap  by  the  reactive  obstacles.  Because  of  the  system’s  symmetry,  cells  arc  not  able  to 
acquire  signal  from  one  another  by  moving  off  the  centerline,  so  cells  are  prevented  from 
reaching  one  another,  in  contrast  to  the  case  with  the  single  obstacle.  This  is  somew'hat 
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surprising:  apparently,  obstacle  boundary  reactivity  in  this  particular  gap  geometry  can 
raise  by  an  order  of  magnitude  the  minimal  gap  size  allowed  for  passage  of  a 
chemotactically  active  cell. 

We  undertook  this  study  as  a first  attempt  to  understand  how  an  imposed  geometrical 
length  scale  (such  as  a pore  size  in  a porous  domain)  together  with  chemoattractant 
transport  with  possible  surface  reactions  influence  chemotactic  cell  behavior.  Clearly,  both 
parameters  are  important,  as  demonstrated  by  figure  5a.  It  remains  to  be  seen  whether 
this  understanding  can  be  translated  into  the  behavior  of  chemotactically  active  cells 
navigating  through  a fully  random  porous  domain  in  three  dimensions,  which  is  our  current 
object  of  study. 

CONCLUSIONS 

We  have  developed  a simulation-based  model  of  the  locomotion  of  cells  which  both 
secrete  and  detect  chemoattractant,  coupled  to  transient  concentration  fields  of  this 
chemoattractant,  for  2D  heterogeneous  domains  with  boundaries  that  consume  (or 
produce)  the  chemoattractant.  This  model  serves  as  a preliminary  step  toward  developing 
a comprehensive  model  of  chemically-stimulated  tissue  ingrowth  in  porous  scaffolds,  and  is 
particularly  suited  to  modeling  angiogenesis.  Clearly,  chemotactic  response  is  not  the  only 
phenomena  governing  tissue  ingrowth;  transport  of  nutrients  and  wastes,  as  well  as,  cell 
adhesion,  spreading,  proliferation  must  eventually  be  taken  into  account.  We  have  simply 
sought  to  understand  better  the  possible  roles  played  by  chemotactive  response  given  that 
(a)  cells  use  diffusing  cytokines  to  signal  one  another,  and  these  signals  are  a crucial 
element  of  tissue  growth,  and  (b)  in  porous  media,  the  chemical  nature  of  the  surfaces  can 
be  engineered.  We  have  found  that  the  presence  of  obstacles  with  reactive  surfaces  greatly 
influences  the  chemotactic  response  of  crawling  cells.  Generally,  reactive  surfaces  repel 
migratory  cells  which  would  normally  respond  by  moving  in  the  direction  of  increasing 
concentration.  We  have  seen  that  moderate  surface  reactivities  are  sufficient  to  prevent 
successful  transmission  of  chemical  signals  through  gaps  much  wider  than  the  minimum 
required  for  cell  transit.  These  results,  while  applicable  only  in  a generic  sense,  hint  at 
possible  mechanisms  for  guiding  the  directionality  and  speed  of  cell  migration  with  fine 
detail,  possibly  offering  guidance  as  to  how  to  engineer  chemically  the  internal  pore 
surfaces  of  a porous  scaffold  to  enhance  the  likelihood  of  successful  tissue  ingrowth. 
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